Lyapunov exponent of the random frequency oscillator: cumulant expansion approach 
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We consider a one-dimensional harmonic oscillator with a random frequency, focusing on both the 
standard and the generalized Lyapunov exponents, A and A* respectively. We discuss the numerical 
difficulties that arise in the numerical calculation of A* in the case of strong intermittency. When 
the frequency corresponds to a Ornstein-Uhlenbeck process, we compute analytically A* by using 
a cumulant expansion including up to the fourth order. Connections with the problem of finding 
an analytical estimate for the largest Lyapunov exponent of a many-body system with smooth 
interactions are discussed. 



I. INTRODUCTION 

The theory of Lyapunov exponents of hard-ball sys- 
tems has a long history. It started with the pioneering 
work of Krylov [l|, Q , was rigorously developed by Sinai 
and collaborators, and completed (to some extent) by 
van Beijercn, Dorfman and co-workers ■ The analyt- 
ical calculation of, e.g., the largest Lyapunov exponent 
of a dilute rigid-sphere gas, is based on the fact that the 
dynamics consists of free rectilinear motions interrupted 
by instantaneous elastic collisions @ ; the expressions so- 
obtained agree quantitatively with the numerical exper- 
iments [III El. 

The case of a dilute gas with finite-range interactions 
can be handled in close analogy with the rigid-sphere 
problem: though the collisions are not trivial any more, 
the dynamics is still ruled by occasional pairwise encoun- 
ters (a. [ill. Il2j . However, when one considers long-range 
interactions (or short-range interactions and high den- 
sities), the theoretical approach must be substantially 
modified. 

In the general case we must deal with the full system 
of coupled differential equations that govern the evolu- 
tion of multidimensional tangent vectors r/(t). Consider 
for instance a gas of N particles in three dimensions de- 
scribed by the Hamiltonian 



3N r> 2 
^ 2m 
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+ V(gi, . . . ,g 3A r), 
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where qi and pi, are conjugate position- momentum co- 
ordinates. Assuming m = 1, tangent vectors evolve ac- 
cording to 
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(dot meaning time derivative), where V is the Hessian 
matrix of the potential V, namely 



d 2 V 
dqidqj 



(3) 



The Hessian depends explicitly on time because it is cal- 
culated along a reference trajectory q(t). Once initial 



conditions x = (qo,Po) an d ?7o have been specified, one 
can find r/(t) from Eq. @. Then the L yap unov exponent 
A is obtained by calculating the limit [13[ 



A = lim ]-ln\ri(t;x ,rio)\ 

t-KX) t 



(4) 



Assuming crgodicity on the energy-shell, A becomes in- 
dependent of initial conditions xq, which can then be 
chosen randomly according to the microcanonical distri- 
bution. There will also be no dependence on initial tan- 
gent vectors, because if ?yo is also chosen randomly, it will 
have a non-zero component along the most expanding di- 
rection. It is this average over Xq and ryo that permits 
to treat equations ^ formally as a system of stochas- 
tic differential equations [lij]. Moreover, if the dynamics 
can be thought of as free motion plus weak interactions, 
then perturbative techniques, like the cumulant expan- 
sion [lJ-Qil, can be invoked. So, the theory attempts to 
calculate the average 



A = lim -(ln\r](t;xo,r]o)\) ■ 

t->oo t 
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However, in practice, it is much simpler to develop an 
estimate for the generalized Lyapunov exponent [l7], 13 
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This is essentially the approach followed by Barnett et 
T9I- 2l| . Pettini et al |22l - [2~ij . and the present authors 
25 27]. In situations of weak intermittency both ex- 
ponents are expected to be close. If one wishes to use 
a theoretically calculated A* as an approximation to A, 
then a numerical check must be done first to verify that 
both exponents coincide. The cumulant expansion to be 
discussed below offers an analytical expression for A via 
the replica trick. Note, however, that the difficulties in- 
volved in such a calculation are much greater than those 
we shall face when dealing with A*. 

Though there are some differences among the formula- 
tions of the three just-mentioned groups, it may be said 
that the main theoretical conclusion extracted from that 
body of work is: if one combines the cumulant expansion 
with some kind of isotropy approximation (which may 
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be fully justified in some cases), the original problem of 
6iV differential equations can be reduced to a system of 
only two equations for a "representative" single degree of 
freedom: 



moments: 



r, 2 J \- K (t) o) U 



(7) 



In this kind of mean- field approximation, the "curvature" 
n(t) is a scalar stochastic process, whose cumulants can 
be related to the (operator) cumulants of the Hessian 
V(i) (see, e.g., 

The comparison of theoretical results obtained with 
the cumulant approach versus numerical simulations has 
met mixed success. The agreement is very good for a 
many-particle system with bounded weak interactions 
[H,|27| and for the Fermi- Past a- Ulam system (24|. How- 
ever, the results for the ld-XY model (24|, for a dense 
one-component plasma (l9l.l28l|. and for a dilute Lennard- 
Jones gas are not so satisfactory. 

The purpose of this paper is to investigate the limits 
of validity of the cumulant approach for the Lyapunov 
exponent of a many particle system. We choose BjS '& 
starting point the simplified mean-field setting ([7]) and 
consider two possibilities for n(t). It has been argued 
that, for typical chaotic many-body systems, K(t) should 
be close to Gaussian white noise; this is the first case 
we shall consider. In the white-noise case the second- 
order expansion for A* is exact, thus this case is ideally 
suited for analyzing the difficulties that appear in the 
numerical calculation of A*. Next, we keep the Gaussian 
and Markov properties but allow for finite correlation 
times, leading to the Ornstein-Uhlenbeck process; in this 
case we calculate the fourth cumulant contribution to A*. 
Though it will not be considered here, we also mention 
the interesting situation of n{t) being a Poisson process, 
which appears to be the appropriate choice for modeling 
the tangent-vector dynamics in a dilute gas with short- 
range interactions. 
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Let us think that, in principle, both parameters a and 
k are stationary stochastic processes. If fluctuations are 
small enough (in a sense that will be discussed later), 
one can obtain the average of the second-moment vector 
using the first terms of the cumulant expansion, which 
works as follows [l4[ • First we split the stochastic matrix 
as an average plus fluctuations: 



B(t) =B +B 1 (t). 
For long times one has: 




(10) 



(11) 



where K is the 3x3 matrix given by the operator cumu- 
lant expansion fl4| 

/>oo 

K = B + / (B 1 (r)e Bor B 1 (0))e- BoT dr + ... . (12) 
Jo 

Dots stand for third and higher cumulants (some explicit 
expressions can be found in fltl]). The exponent A* is 
related to the eigenvalue of K that has the largest real 
part: 



A* = — max -ft {fci, fc2j ^3} 



(13) 



with ki the eigenvalues of K. 



III. GAUSSIAN WHITE NOISE 



II. CUMULANT EXPANSION FOR THE KUBO 
OSCILLATOR 

Formally, Eq. ([7J describes a harmonic oscillator with a 
random frequency u) such that cu 2 = k (Kubo oscillator). 
It is worth generalizing this model a bit to account for the 
possibility of damping, i.e., we shall consider an oscillator 
described by the dynamical equation 

q + aq + nq = . (8) 

Setting a = 0, q = 771, p = q = 772, recovers ([7]). 

Some analytical results for the Lyapunov exponent of 
the Kubo oscillator ([8]) can be found in the literature 
(see, e.g., [3(1 HJ). Here we shall restrict ourselves to 
the analytical calculation of the generalized exponent A*. 
For this purpose we must consider the dynamics of second 



When the entries of the fluctuation matrix Bi are 
Gaussian white noise (and only in this case [l6| ) the cu- 
mulant expansion stops at the second order, i.e., Eq. (TT2")) 
is exact (without the ellipsis). This is the case we con- 
sider now. 



A. Random frequency 

Let us first study the situation where the damping a 
is a constant and 

K(t) = Ko +at), (14) 

where £(t) is a zero-mean Gaussian white noise. Its cor- 
relation function reads 

(f(t)f(0)=A*(i-0- (15) 
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With these definitions one has 



B= -2a -2k +£(*) 2. (16) 
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After substitution into Eq. (fT2")l we readily obtain 



K = 



The generalized exponent A* can now be calculated from 
Eq. (fl~3| . A closed expression for the standard Lyapunov 
exponent can be found in the literature [30j ]. As an ex- 
ample, Fig. [T] displays analytical results for both expo- 
nents. We also show the outcomes of numerical simula- 
tions. Given that the theoretical expressions are exact, 
Fig. [T] constitutes a test for our numerical calculations. 
Numerical details, including a discussion about the diffi- 
culties found in the calculation of A*, will be presented 
in Sec. HVl 




([HI) will be taken in Stratonovich sense. Therefore, 
matrix B in Eq. © can be decomposed as 
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Hence, substitution into (fT2j) yields 
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(19) 



(20) 



From the eigenvalues of K we obtain A* following 
Eq. (1131) . A theoretical expression for A can be found in 
Ref. [3l|. Fig. [5] exhibits numerical and analytical results 
for both exponents, as a function of the noise intensity. 




FIG. 2: Harmonic oscillator with random damping. Shown 
are the Lyapunov exponents A and A* as a function of the 
noise strength A. Solid lines correspond to the theoretical 
expressions given by (|13I20[) and Ref. [3l| . for A* and A, re- 
spectively. Symbols are the results of numerical simulations 
(averaged over 10 4 trajectories). We chose k — 1 and ao = 1. 



FIG. 1: Harmonic oscillator with random frequency. Shown 
are the Lyapunov exponents A and A* as a function of the 
noise strength A. Solid lines correspond to the theoretical 
expressions given by inserting (Ti7|) into (Tl3jl , and in Ref. [33] 
for A* and A, respectively. Symbols are the results of numer- 
ical simulations (averaged over 10 4 trajectories). We chose 
a = and kq = 1. 



B. Random damping 

Now we consider an harmonic oscillator with constant 
frequency but in an environment with fluctuating damp- 
ing coefficient 

£*(«) = a +£(*) > (18) 

where is also in this case a zero-mean Gaussian white 
noise. The corresponding stochastic differential equation 



Note that in the cases considered above A and A* do 
not coincide. We have checked that the difference be- 
tween them (which is a quantifier of the degree of inter- 
mittency of the dynamics) may be controlled by suitable 
choice of the parameters of the oscillator. We preferred to 
consider intermittent cases, because it is in these regimes 
that the numerical difficulties arise, as we discuss in the 
following section. We remark that there are situations of 
interest where both exponents practically coincide, e.g., 
for a dilute Lennard- Jones gas |29|. In such cases a the- 
ory capable of estimating A* will also produce a good 
estimate for the standard Lyapunov exponent A. 

IV. NUMERICAL TREATMENT 

Numerical simulations were performed by means of the 
Euler algorithm with time step dt = 10 -3 . For each tra- 
jectory we computed the norm |?7(£)| 2 = q 2 +p 2 as a func- 
tion of time t. The Lyapunov exponent is approximated 
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by the average over initial conditions of the finite-time 
exponents: 

A -<ln |T7(*;ar ,T7o)|> = (X(t;x ,r]o)) , (21) 

where t is large enough to guarantee the convergence of 
the average to the desired precision. In order to obtain 
the asymptotic value of the generalized exponent ([6]), in 
principle, one must calculate the squared-norm averaged 
over several realizations at a given large time. However, 
we must keep in mind that such an average is dominated 
by the extreme positive values of the local exponent X(t). 
Hence, direct averaging over n(i)| 2 may yield spurious 
results whenever the variance of X(t) fails to vanish with 
time fast enough. To avoid this problem, instead of the 
straightforward averaging, we preferred to estimate the 
local generalized exponent from the cumulants of the dis- 
tribution of X(t): 



X*(t) 



ln(\v(t)\ 2 ) _ ln(e 2A (*)*) 



2t 



2t 



E 

n>l 



(22) 

where the nth-order cumulants of the distribution 

of local exponents X(t). Fig. [3] illustrates, for the white- 
noise random frequency oscillator, the (X(t)) as a func- 
tion of time (first-order truncation of (f2"2"j) ) , as well as the 
expansion (|22[) truncated at the second and third orders. 
For comparison, also plotted is the crude estimate (J5]). 
Clearly, the expansion ([22j) has to be considered in order 
to properly estimate A*. In Figs. [T] and [21 A* was numeri- 
cally computed from the third-order truncation, because 
the next (noisier) terms do not contribute significantly. 




40 60 

t 



FIG. 3: Numerical difficulties in the calculation of A*. We 
plot the finite-time exponent A(t) = «i(i) as a function of 
time (red line) for the random frequency oscillator with a — 
and A = 50. Averages were computed over 10 5 realizations. 
Also plotted are the corrections arising from the second and 
third cumulants of the distribution of finite-time Lyapunov 
exponents Eq. (|22[) (blue). For comparison we also show 
the straightforward average ln(|rj(t; xo, rjo)\ 2 )/2t (light gray). 
Dashed lines correspond to the theoretical asymptotic values. 



V. CORRELATED NOISE 

For white noise fluctuations, either in the frequency or 
in the damping, we have verified in Sec. 3 (see Figs. [T] 
and [2]) , that our theory for A* is in agreement with nu- 
merical results, provided the later are obtained by means 
of the procedure described in the preceding section. The 
analysis in Sec. 3 also allows to quantify the discrepancy 
between A and A*, which typically increases with increas- 
ing amplitude of the fluctuations. 

Now we shall analyze the effect of introducing noise 
correlations. We consider again the case of a random 
frequency, as in Eq. ()14|) . but now the noise is a zero- 
mean Ornstein-Ulhenbeck process, i.e., with correlation 
function 



mat')) 



A 
27 



exp(-|t - t'\/r) = cr 2 exp(-|t-i'|/ T ) 



(23) 

For simplicity we set a = and kq = 0. By inserting 
Eq. (|16|) into Eq. (fT2|) . the second-cumulant matrix K^ 2 ' 
becomes 



K (2) 




(24) 



Notice that in the limit r — > the white-noise case is 
recovered. 

In the presence of correlations the second-order trun- 
cation of the cumulant expansion (fl2|) is not exact. In 
order to improve the theory one must calculate higher 
cumulants. For the present case the third cumulant is 
null. Explicit expressions for the fourth cumulant were 
given by Fox [l6[ and by Breuer et al (32|. So, the fourth 
cumulant can be calculated without great effort (with 
the aid of algebraic manipulation programs) . The fourth 
order approximation to K reads 



KW = K< 2 > + ±AV 




(25) 



The comparison between the theoretical results for A* 
(with the second (blue) and fourth (dark blue) order cor- 
rections) and numerical outcomes is shown in Fig. 01 No- 
tice that in numerical estimates the fourth-order correc- 
tion is very small in comparison with the third-order one, 
suggesting that the cumulant expansion is rapidly con- 
verging. 



A. Kubo number 

The perturbation parameter controlling the conver- 
gence of the cumulant expansion is the so-called Kubo 
number e. General considerations led van Kampen [Tij to 
conclude that the Kubo number is the product of the am- 
plitude of the fluctuations and the correlation time, that 
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FIG. 4: Harmonic oscillator with correlated random fre- 
quency. We show A and A* as a function of the noise ampli- 
tude A. Solid lines correspond to the theoretical results ob- 
tained from (|13I24|) (blue) and (|13l25j) (dark blue) for A* and 
to the approximate expression (using the decoupling ansatz) 
following Ref. [3(| for A. Symbols are the results of numerical 
simulations (averaged over 10 5 trajectories), corresponding to 
the second (circles) , third (squares) and fourth (triangles) or- 
der corrections of (|22[) . Parameters are a = 0, kq = and 

T = 1. 

is err. However, in the present case it is clear that such 
a combination is not adimensional. The correct Kubo 
number is instead 




This can be checked explicitly from the second and fourth 
cumulants above. Consider, for instance, the element 
K21, which dominates the Lyapunov exponent for small 
correlation times: 

k 21 = a + H a v + ... = a(\ + H At 3 + ..^ . 

(27) 



In the white-noise limit, i.e., r — > with A fixed, the 
Kubo number tends to zero -as it should be. 



VI. FINAL REMARKS 

We have taken the first step towards the application 
of the cumulant expansion to calculate the largest Lya- 
punov exponent of a dilute gas. 

The case of white-noise fluctuations (either in the fre- 
quency or in the damping) was considered first. This 
study was very useful to understand the difficulties be- 
hind the numerical calculation of the generalized expo- 
nent A*. It was verified that A* can be obtained with 
a satisfactory precision by using the cumulant expansion 
for the distribution of the finite-time Lyapunov exponent. 

We also analyzed briefly the case of correlated noise. 
For the Ornstein-Uhlcnbeck noise we were able to ob- 
tain the fourth cumulant contribution to the analytical 
A*, which showed an improvement with respect to the 
second order truncation, when compared with numerical 
outcomes. Moreover, we showed that the correct pertur- 
bative parameter for the present problem is the product 
err 2 , and not err, as a literal reading of van Kampcn's 
discussion [l4| would suggest. 

It is expected that the present results will be helpful 
for the correct application of the cumulant approach in 
higher dimensionality systems, as well as for the numer- 
ical checking of its validity. 
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